{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `fit_residuals.ipynb`\n",
    "Notebook by Tim Bartholomaus, June 18, 2019\n",
    "\n",
    "Goal of notebook is to plot residuals between a reference profile/surface as a function of space and time, and then fit a smooth, interpolated surface through these residuals.\n",
    "\n",
    "Residuals are expected in a dictionary, `res` with keys `\"x\"`, `\"t\"`, and `\"off\"`, in which `x` is position along the profile, `t` is the relative time of the residual measurement, and `off` is the measured offset between the ICESat-2 measurement and the reference surface elevation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "from scipy import interpolate\n",
    "\n",
    "import sys\n",
    "sys.path.insert(0, '~/gridding/notebook')\n",
    "# import utils\n",
    "\n",
    "# print(sys.path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['~/gridding/notebook', '/srv/conda/lib/python36.zip', '/srv/conda/lib/python3.6', '/srv/conda/lib/python3.6/lib-dynload', '', '/srv/conda/lib/python3.6/site-packages', '/srv/conda/lib/python3.6/site-packages/IPython/extensions', '/home/jovyan/.ipython']\n"
     ]
    }
   ],
   "source": [
    "print(sys.path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create synthetic residuals that are a function of length along the profile `x`, and time `t`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# n_residuals = 30\n",
    "# max_length = 50*1000 # m\n",
    "# max_time = 100 # d\n",
    "# max_off = 100 # m\n",
    "\n",
    "\n",
    "# res = dict()\n",
    "# res['x'] = max_length * np.random.rand(n_residuals) # m\n",
    "# res['t'] = max_time * np.random.rand(n_residuals) # d\n",
    "\n",
    "# noise = np.random.randn(n_residuals) * 20 # m\n",
    "\n",
    "# off_scalar = max_off / (max_length*max_time)\n",
    "# res['off'] = off_scalar * res['x']*res['t'] + noise # create the synthetic residuals, with noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read in measured residuals that are a function of length along the profile `x`, and time `t`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "atm_date = '2014-04-29'\n",
    "# atm_date = '2018-04-18'\n",
    "\n",
    "data_prod_path = '../data_prod/'\n",
    "# intersections_filename = 'Intersections_ATM20180418.csv'\n",
    "# intersections_filename = 'Intersections_ATM'+atm_date[0:4] + atm_date[5:7] + atm_date[8:10] + '.csv'\n",
    "intersections_filename = 'residuals.csv'\n",
    "res = pd.read_csv(data_prod_path + intersections_filename, parse_dates=[5]) \n",
    "\n",
    "res.head()\n",
    "\n",
    "res.dtypes\n",
    "\n",
    "res['off'] = res['residuals']\n",
    "# res['off'] = res['z_ATL06'] - res['ATM_elev']\n",
    "res['t'] = res['t_ATL06'] - pd.to_datetime(atm_date)\n",
    "res['t_days'] = res['t'].astype('timedelta64[D]')\n",
    "res['t_flt'] = res['t_days'].astype('float64')\n",
    "res.rename(columns={'dist_along': \"x\"}, inplace=True)\n",
    "\n",
    "\n",
    "# res.head()\n",
    "# # res.columns\n",
    "\n",
    "# pd.to_datetime(atm_date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res['t_ATL06'] - pd.to_datetime(atm_date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# n_residuals = 30\n",
    "# max_length = 50*1000 # m\n",
    "# max_time = 100 # d\n",
    "# max_off = 100 # m\n",
    "\n",
    "\n",
    "# res = dict()\n",
    "# res['x'] = max_length * np.random.rand(n_residuals) # m\n",
    "# res['t'] = max_time * np.random.rand(n_residuals) # d\n",
    "\n",
    "# noise = np.random.randn(n_residuals) * 20 # m\n",
    "\n",
    "# off_scalar = max_off / (max_length*max_time)\n",
    "# res['off'] = off_scalar * res['x']*res['t'] + noise # create the synthetic residuals, with noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Spatio-temporal (2-D) interpolation\n",
    "Right now, the radial basis functions seem to perform best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     1633 days 15:53:52\n",
       "1     1633 days 15:53:52\n",
       "2     1633 days 15:53:52\n",
       "3     1633 days 15:53:52\n",
       "4     1636 days 05:21:45\n",
       "5     1640 days 05:13:27\n",
       "6     1640 days 05:13:27\n",
       "7     1640 days 05:13:27\n",
       "8     1640 days 05:13:27\n",
       "9     1640 days 05:13:27\n",
       "10    1641 days 15:37:16\n",
       "11    1641 days 15:37:16\n",
       "12    1641 days 15:37:16\n",
       "13    1641 days 15:37:16\n",
       "14    1641 days 15:37:16\n",
       "15    1641 days 15:37:16\n",
       "16    1657 days 04:14:44\n",
       "17    1657 days 04:14:44\n",
       "18    1657 days 04:14:44\n",
       "19    1657 days 04:14:44\n",
       "20    1657 days 04:14:44\n",
       "21    1657 days 04:14:44\n",
       "22    1661 days 04:06:22\n",
       "23    1661 days 04:06:22\n",
       "24    1661 days 04:06:22\n",
       "25    1661 days 04:06:22\n",
       "26    1661 days 04:06:22\n",
       "27    1661 days 04:06:22\n",
       "28    1665 days 03:57:57\n",
       "29    1665 days 03:57:57\n",
       "             ...        \n",
       "153   1747 days 23:54:28\n",
       "154   1747 days 23:54:28\n",
       "155   1747 days 23:54:28\n",
       "156   1747 days 23:54:28\n",
       "157   1751 days 23:46:06\n",
       "158   1751 days 23:46:06\n",
       "159   1751 days 23:46:06\n",
       "160   1751 days 23:46:06\n",
       "161   1751 days 23:46:06\n",
       "162   1751 days 23:46:06\n",
       "163   1753 days 10:09:55\n",
       "164   1753 days 10:09:55\n",
       "165   1755 days 23:37:46\n",
       "166   1755 days 23:37:46\n",
       "167   1755 days 23:37:46\n",
       "168   1755 days 23:37:46\n",
       "169   1755 days 23:37:46\n",
       "170   1755 days 23:37:46\n",
       "171   1757 days 10:01:34\n",
       "172   1757 days 10:01:34\n",
       "173   1757 days 10:01:34\n",
       "174   1757 days 10:01:34\n",
       "175   1757 days 10:01:34\n",
       "176   1757 days 10:01:34\n",
       "177   1761 days 09:53:16\n",
       "178   1761 days 09:53:16\n",
       "179   1761 days 09:53:16\n",
       "180   1761 days 09:53:16\n",
       "181   1761 days 09:53:16\n",
       "182   1761 days 09:53:16\n",
       "Name: t, Length: 183, dtype: timedelta64[ns]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res['t']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0      122581.170689\n",
       "1      122647.058967\n",
       "2      125810.129706\n",
       "3      125875.304455\n",
       "4       77393.293766\n",
       "5       67655.902876\n",
       "6       67687.240097\n",
       "7       64484.432996\n",
       "8       61379.685494\n",
       "9       61315.684441\n",
       "10      89377.265064\n",
       "11      89479.464606\n",
       "12      92780.627483\n",
       "13      92848.522698\n",
       "14      96176.153026\n",
       "15      96241.609114\n",
       "16     121954.028971\n",
       "17     121854.851720\n",
       "18     118653.494370\n",
       "19     118588.415511\n",
       "20     115361.066547\n",
       "21     115295.191381\n",
       "22     105596.029174\n",
       "23     105001.339697\n",
       "24     102145.248209\n",
       "25     102080.404536\n",
       "26      98960.114430\n",
       "27      98894.603155\n",
       "28      89241.058144\n",
       "29      89172.935318\n",
       "           ...      \n",
       "153    121425.520431\n",
       "154    121359.531822\n",
       "155    118164.360812\n",
       "156    118098.943857\n",
       "157    108153.357245\n",
       "158    108088.013034\n",
       "159    104835.607040\n",
       "160    104769.182503\n",
       "161    101658.707466\n",
       "162    101593.796438\n",
       "163    125289.795232\n",
       "164    125354.749926\n",
       "165     91929.771598\n",
       "166     91861.692774\n",
       "167     88729.738750\n",
       "168     88661.510375\n",
       "169     85556.001064\n",
       "170     85454.321057\n",
       "171    109036.355948\n",
       "172    109134.723785\n",
       "173    112275.636373\n",
       "174    112341.684666\n",
       "175    115492.850024\n",
       "176    115591.706990\n",
       "177     92270.364511\n",
       "178     92338.479124\n",
       "179     95651.996664\n",
       "180     95717.578915\n",
       "181     98960.114430\n",
       "182     99025.590148\n",
       "Name: x, Length: 183, dtype: float64"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res['x']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Create a mesh of the space and time variables\n",
    "# grid_x, grid_t = np.mgrid[0:max_length:100j, \n",
    "#                           0:max_time:100j] # create a grid from 0 to max_length and 0 to\n",
    "#                                            #  max_time, with the j number of steps\n",
    "# Create a mesh of the space and time variables\n",
    "grid_x, grid_t = np.mgrid[50000:130000:100j, \n",
    "                          1620:1780:100j] # create a grid from 0 to max_length and 0 to\n",
    "                                           #  max_time, with the j number of steps\n",
    "    \n",
    "# # # This griddata is too strict\n",
    "grid_off = interpolate.griddata((res['x'].values, res['t_flt'].values), res['off'].values, \n",
    "                                (grid_x, grid_t), method='linear')\n",
    "\n",
    "# tck = interpolate.bisplrep(res['x'].values, res['t'].values, res['off'].values, s=2)\n",
    "# grid_off = interpolate.bisplev(grid_x[:,0], grid_t[0,:], tck)\n",
    "\n",
    "# f = interpolate.interp2d(res['x'].values, res['t_flt'].values, res['off'].values, kind='linear')\n",
    "# grid_off = f(grid_x[:,0], grid_t[0,:])\n",
    "\n",
    "# # Radial basis functions\n",
    "# x_stan = res['x']/150#max_length\n",
    "# t_stan = res['t']/1800#max_time\n",
    "# rbf = interpolate.Rbf(x_stan, t_stan, res['off'], epsilon=1)\n",
    "#         # epsilon defines the stiffness of the fitting routine. High values are stiff surfaces\n",
    "# grid_x_stan = grid_x/150#max_length\n",
    "# grid_t_stan = grid_t/1800#max_time\n",
    "# grid_off = rbf(grid_x_stan, grid_t_stan)\n",
    "\n",
    "\n",
    "# # Kriging\n",
    "# # Function does not like NaN's, so lets remove them\n",
    "# # xp, yp, zp = x[~np.isnan(z)], y[~np.isnan(z)], z[~np.isnan(z)]\n",
    "# Run kriging interpolator, with 1 m RMSE noise added to the diagonal (np.ones(xp.shape)*1)\n",
    "# grid_off_lsc, e_lsc = utils.lscip(res['x'], res['t'], res['off'], np.ones(res['x'].shape)*1, grid_x, grid_t, n=10, d=1e3, a=25e3)[0:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan, nan, nan, ..., nan, nan, nan],\n",
       "       [nan, nan, nan, ..., nan, nan, nan],\n",
       "       [nan, nan, nan, ..., nan, nan, nan],\n",
       "       ...,\n",
       "       [nan, nan, nan, ..., nan, nan, nan],\n",
       "       [nan, nan, nan, ..., nan, nan, nan],\n",
       "       [nan, nan, nan, ..., nan, nan, nan]])"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_off"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3c2cf90c6dc84837a5f5b5e470a93df3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "plt.figure(figsize=(10,6))\n",
    "# plt.contourf(grid_x/1000, grid_t, (grid_off))\n",
    "plt.scatter(res['x'].values/1000, res['t_days'].values, s = 20, c=res['off'].values,\n",
    "            linewidth=.5, edgecolor='k', vmin = -50, vmax=10, cmap='inferno')\n",
    "# res.plot.scatter(x='x', y='t', c='off')#, linewidth=.5, edgecolor='k')\n",
    "plt.xlabel('Profile Length (km)')\n",
    "plt.ylabel('Time (yrs)')\n",
    "plt.xlim(50,130)\n",
    "plt.ylim(1620,1780)\n",
    "# plt.colormap('RdBu')\n",
    "plt.colorbar(label='Deviation from reference profile (m)');\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 50.        ,  50.        ,  50.        , ...,  50.        ,\n",
       "         50.        ,  50.        ],\n",
       "       [ 50.80808081,  50.80808081,  50.80808081, ...,  50.80808081,\n",
       "         50.80808081,  50.80808081],\n",
       "       [ 51.61616162,  51.61616162,  51.61616162, ...,  51.61616162,\n",
       "         51.61616162,  51.61616162],\n",
       "       ...,\n",
       "       [128.38383838, 128.38383838, 128.38383838, ..., 128.38383838,\n",
       "        128.38383838, 128.38383838],\n",
       "       [129.19191919, 129.19191919, 129.19191919, ..., 129.19191919,\n",
       "        129.19191919, 129.19191919],\n",
       "       [130.        , 130.        , 130.        , ..., 130.        ,\n",
       "        130.        , 130.        ]])"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b18d1b763d4e40338afecf5c2a8977f0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# fig, ax = plt.subplots(3)\n",
    "\n",
    "# plt.figure()\n",
    "plt.figure(figsize=(10,6))\n",
    "\n",
    "plt.contourf(grid_x/1000, grid_t, grid_off, cmap='inferno')\n",
    "# plt.scatter(res['x']/1000, res['t'], c=res['off'], linewidth=.5, edgecolor='k')\n",
    "plt.scatter(res['x'].values/1000, res['t_days'].values, s = 20, c=res['off'].values,\n",
    "            linewidth=.5, edgecolor='k', vmin = -50, vmax=10, cmap='inferno')\n",
    "plt.xlabel('Profile Length (km)')\n",
    "plt.ylabel('Time (days)')\n",
    "plt.colorbar(label='Deviation from reference profile (m)');\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1620.         1621.61616162 1623.23232323 ... 1776.76767677\n",
      "  1778.38383838 1780.        ]\n",
      " [1620.         1621.61616162 1623.23232323 ... 1776.76767677\n",
      "  1778.38383838 1780.        ]\n",
      " [1620.         1621.61616162 1623.23232323 ... 1776.76767677\n",
      "  1778.38383838 1780.        ]\n",
      " ...\n",
      " [1620.         1621.61616162 1623.23232323 ... 1776.76767677\n",
      "  1778.38383838 1780.        ]\n",
      " [1620.         1621.61616162 1623.23232323 ... 1776.76767677\n",
      "  1778.38383838 1780.        ]\n",
      " [1620.         1621.61616162 1623.23232323 ... 1776.76767677\n",
      "  1778.38383838 1780.        ]] [1633. 1633. 1633. 1633. 1636. 1640. 1640. 1640. 1640. 1640. 1641. 1641.\n",
      " 1641. 1641. 1641. 1641. 1657. 1657. 1657. 1657. 1657. 1657. 1661. 1661.\n",
      " 1661. 1661. 1661. 1661. 1665. 1665. 1665. 1665. 1665. 1665. 1666. 1666.\n",
      " 1666. 1666. 1666. 1666. 1669. 1669. 1669. 1669. 1669. 1669. 1670. 1670.\n",
      " 1670. 1670. 1670. 1670. 1674. 1674. 1674. 1674. 1674. 1674. 1686. 1686.\n",
      " 1686. 1686. 1686. 1686. 1690. 1690. 1690. 1690. 1690. 1690. 1694. 1694.\n",
      " 1694. 1694. 1694. 1694. 1695. 1695. 1695. 1695. 1695. 1695. 1698. 1698.\n",
      " 1698. 1698. 1698. 1698. 1699. 1699. 1699. 1699. 1699. 1699. 1702. 1702.\n",
      " 1702. 1702. 1703. 1703. 1703. 1703. 1703. 1703. 1707. 1707. 1707. 1707.\n",
      " 1707. 1707. 1719. 1719. 1719. 1719. 1719. 1723. 1723. 1723. 1723. 1723.\n",
      " 1723. 1724. 1724. 1724. 1724. 1724. 1724. 1727. 1727. 1727. 1727. 1727.\n",
      " 1727. 1731. 1731. 1731. 1731. 1731. 1731. 1732. 1732. 1732. 1732. 1732.\n",
      " 1732. 1736. 1736. 1736. 1736. 1736. 1736. 1747. 1747. 1747. 1747. 1747.\n",
      " 1747. 1751. 1751. 1751. 1751. 1751. 1751. 1753. 1753. 1755. 1755. 1755.\n",
      " 1755. 1755. 1755. 1757. 1757. 1757. 1757. 1757. 1757. 1761. 1761. 1761.\n",
      " 1761. 1761. 1761.]\n"
     ]
    }
   ],
   "source": [
    "print(grid_t, res['t_days'].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid_t[:,50]\n",
    "np.savetxt('residuals_at_1700d.csv', grid_off[:,50]) \n",
    "np.savetxt('along_dist_at_1700d.csv', grid_x[:,50]) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Timestamp('2018-12-24 06:05:01.000030')"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print(pd.to_datetime(atm_date) + pd.Timedelta('0')#1700:00:00.000') )\n",
    "pd.to_datetime(atm_date) + pd.to_timedelta('1700 days 06:05:01.00003')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
